%-------------------------------------------------------------------------;
% Table 6: Predicting the Term Structure in Market Model Alphas;
%-------------------------------------------------------------------------;

clear;
clc;
close all;

%-------------------;
% Import data;
%-------------------;

filename='Main_data_1996_2022';
[data,txt,raw]=xlsread(filename);

dates=txt(2:end,1);

%Index data;
SP_500=data(:,1);         %SP 500 index price level;
Div_m=data(:,2);            %Monthly index dividends;
Div_a=data(:,3);            %Annual index dividends;
SP_ret=data(2:end,4);       %Monthly return on SP 500 index;

%Dividend strip data;
Div_strip_12=data(:,5);     %Price of 12-month dividend strip;
Div_strip_ret=data(2:end,6);%Monthly return on dividend strip strategy;

%Indicated dividends;
Ind_div=data(1:end,7);  

%Bond returns;
Bond_2y=data(2:end,8);      %Monthly return on 2-year Treasury bond;
Bond_10y=data(2:end,9);     %Monthly return on 10-year Treasury bond;

%Fama-French data;
Mkt_Rf=data(2:end,10);      %Market factor;
Rf=data(2:end,11);          %Risk-free rate;

%Define dividend-to-price ratios;
DP_sp=Div_a./SP_500;
DP_strip=Div_a./Div_strip_12;

%Define indicated dividend growthh;
Ind_dg=(Ind_div./Div_a)-1;     

%Take logarithms of all the main variables;
ret_sp=log(SP_ret+1);
ret_strip=log(Div_strip_ret+1);

bond_2y=log(Bond_2y+1);  
bond_10y=log(Bond_10y+1);

mkt_rf=log(Mkt_Rf+1);
rf=log(Rf+1);

dp_sp=log(DP_sp);
dp_strip=log(DP_strip);

ind_dg=log(Ind_dg+1);

%--------------------------------------------------;
%Rolling window Sharpe ratios and alphas;
%--------------------------------------------------;

%Calculate rho;
rho=(exp(-mean(dp_sp)))/(1+exp(-mean(dp_sp)));

%Take average for dp_strip and start the sample with observation no. 3;
dp_strip_n=filter(ones(3,1),1,dp_strip)./3;

dp_sp=dp_sp(3:end);
dp_strip=dp_strip_n(3:end);

ret_sp=ret_sp(3:end);
ret_strip=ret_strip(3:end);
rf=rf(3:end);

%To store results;
Tab_6=zeros(10,36);

%Calculate excess returns;
ret_sp=ret_sp-rf;
ret_strip=ret_strip-rf;

for h=1:36

    %To store results;
    ret_sp_h=zeros(size(ret_sp,1),1);
    ret_strip_h=zeros(size(ret_strip,1),1);
    std_sp_h=zeros(size(ret_sp,1),1);
    std_strip_h=zeros(size(ret_strip,1),1);
    alpha_h=zeros(size(ret_sp,1),1);

    for i=h:size(ret_strip,1)
        ret_sp_h(i)=sum(ret_sp(i-h+1:i));
        std_strip_h(i)=std(ret_strip(i-h+1:i));
        ret_strip_h(i)=sum(ret_strip(i-h+1:i));
        std_sp_h(i)=std(ret_sp(i-h+1:i));

        y=ret_strip_h;
        x=[ones(size(ret_strip_h,1),1) ret_sp_h];
        [b,~,~,~,~,~]=olsnw(y,x,1,h+6);
        alpha_h(i)=b(1,1);
    end

ret_strip_h=ret_strip_h(h:end);
ret_sp_h=ret_sp_h(h:end);
std_strip_h=std_strip_h(h:end);
std_sp_h=std_sp_h(h:end);
alpha_h=alpha_h(h:end);

SR_strip_h=ret_strip_h./std_strip_h;
SR_sp_h=ret_sp_h./std_sp_h;

%Number of lags for Newey-West for overlapping observations;
lags_overlap=h;

%Number of lags for Newey-West for nonoverlapping observations;
lags_nonoverlap=3;

%Adjust the length of the data to given h; 
dp_sp_h=dp_sp(1:(size(ret_sp_h,1)));
dp_strip_h=dp_strip(1:(size(ret_sp_h,1)));

%For non-overlapping regressions;
t=h;

%------------------------------;
% Table 6;
%------------------------------;

%Difference in alphas;
y=-alpha_h; 

%------------------------------;
%Dp_sp;
%------------------------------;

x=[ones(size(ret_sp_h,1),1) dp_sp_h];
[b,~,tb,R2,~,~]=olsnw(y,x,1,lags_overlap);

Tab_6(1:2,h)=[b(2,1);tb(2,1)];
Tab_6(4,h)=R2;

%To store results;
b_non=zeros(2,h);
tb_non=zeros(2,h);

for i=1:h
    y_non=y(i:h:end);
    x_non=x(i:h:end,:);
    [b,~,tb,~,~,~]=olsnw(y_non,x_non,1,lags_nonoverlap);
    b_non(:,i)=b;
    tb_non(:,i)=tb;
end

t_non=mean(tb_non,2);
Tab_6(3,h)=t_non(2,1);

%------------------------------;
%Scaled dp_mkt - dp_strip;
%------------------------------;

%The AR(1) coefficient for the persistence of expected returns is obtained
%from auxiliary regressions of market returns on dividend-price ratio at
%different horizons (details in the paper).;

AR1=0.851;
A=(1-rho*(AR1^(h/12)));

x=[ones(size(ret_sp_h,1),1) A*dp_sp_h-dp_strip_h];
[b,~,tb,R2,~,error]=olsnw(y,x,1,lags_overlap);

Tab_6(6:7,h)=[b(2,1);tb(2,1)];
Tab_6(9,h)=R2;
Tab_6(10,h)=length(error);

%To store results;
b_non=zeros(2,h);
tb_non=zeros(2,h);

for i=1:h
    y_non=y(i:h:end);
    x_non=x(i:h:end,:);
    [b,~,tb,~,~,~]=olsnw(y_non,x_non,1,lags_nonoverlap);
    b_non(:,i)=b;
    tb_non(:,i)=tb;
end

t_non=mean(tb_non,2);
Tab_6(8,h)=t_non(2,1);
end

Tab_6=Tab_6(1:end,[12 18 24 30 36]);

disp('Table 6')
disp(Tab_6)

